Plot Data
library(ggplot2)
# raw data
ggplot(dataset, aes(x=H2O2, y=Counts)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=siRNA)) +
geom_point(aes(colour=siRNA, shape=Experiment), size=2) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
scale_shape_manual(values=15:20) +
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)")+
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)")+
scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour=siRNA), size=2) +
geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
scale_color_manual(values=c("#000000","#FF0000"))

pdf("FigureS3_H2O2.pdf", width = 5, height = 4)
ggplot(dataset, aes(x=H2O2, y=NormCounts2)) +
theme_bw() +
theme(panel.grid=element_blank(), text = element_text(size=14)) +
geom_point(aes(colour = siRNA, shape = genotype), size=1.75) +
geom_smooth(method=lm, formula = y ~ poly(x,2), se=TRUE,
aes(group = GSID,colour = siRNA, linetype = genotype), fill='#DDDDDD', size=0.5) +
#facet_grid(. ~ genotype) +
xlab(label = "H2O2 (log10 µM)") +
ylab(label = "Normalized Counts") +
scale_color_manual(values=c("#000000","#FF0000")) +
guides(linetype = guide_legend(override.aes= list(color = "#555555")))
dev.off()
## quartz_off_screen
## 2
Models
library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)
Linear formula
fit1 <- lm(Counts ~ Experiment + H2O2*siRNA*genotype, data = dataset)
print(summary(fit1))
##
## Call:
## lm(formula = Counts ~ Experiment + H2O2 * siRNA * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -697.48 -230.97 3.66 207.16 711.48
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2008.20 197.38 10.174 2.11e-12 ***
## Experimentexp2 148.69 127.73 1.164 0.2516
## Experimentexp3 -9.75 127.73 -0.076 0.9396
## H2O2 -1018.50 155.34 -6.557 9.85e-08 ***
## siRNAsiPARP1 -117.27 258.93 -0.453 0.6532
## genotypeALC1KO 361.79 258.93 1.397 0.1704
## H2O2:siRNAsiPARP1 -116.96 219.68 -0.532 0.5975
## H2O2:genotypeALC1KO -410.10 219.68 -1.867 0.0697 .
## siRNAsiPARP1:genotypeALC1KO 423.11 366.19 1.155 0.2551
## H2O2:siRNAsiPARP1:genotypeALC1KO -91.64 310.68 -0.295 0.7696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 361.3 on 38 degrees of freedom
## Multiple R-squared: 0.8871, Adjusted R-squared: 0.8603
## F-statistic: 33.17 on 9 and 38 DF, p-value: 2.56e-15
cat("AIC: ", AIC(fit1))
## AIC: 712.4085
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit2))
##
## Call:
## lm(formula = NormCounts ~ H2O2 * siRNA * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66792 -0.08492 0.07217 0.16452 0.48643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9297 0.1435 13.447 < 2e-16 ***
## H2O2 -0.9597 0.1218 -7.882 1.14e-09 ***
## siRNAsiPARP1 0.3793 0.2029 1.869 0.0690 .
## genotypeALC1KO 0.4093 0.2029 2.017 0.0505 .
## H2O2:siRNAsiPARP1 -0.3915 0.1722 -2.274 0.0284 *
## H2O2:genotypeALC1KO -0.4225 0.1722 -2.454 0.0186 *
## siRNAsiPARP1:genotypeALC1KO -0.3275 0.2870 -1.141 0.2606
## H2O2:siRNAsiPARP1:genotypeALC1KO 0.3381 0.2435 1.388 0.1727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2832 on 40 degrees of freedom
## Multiple R-squared: 0.9189, Adjusted R-squared: 0.9047
## F-statistic: 64.75 on 7 and 40 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit2))
## AIC: 24.33749
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit3))
##
## Call:
## lm(formula = NormCounts2 ~ H2O2 * siRNA * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25608 -0.04318 0.02865 0.07163 0.18649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.98133 0.05621 17.457 < 2e-16 ***
## H2O2 -0.48804 0.04769 -10.233 9.9e-13 ***
## siRNAsiPARP1 -0.09607 0.07950 -1.208 0.234
## genotypeALC1KO -0.08123 0.07950 -1.022 0.313
## H2O2:siRNAsiPARP1 -0.03001 0.06745 -0.445 0.659
## H2O2:genotypeALC1KO -0.04386 0.06745 -0.650 0.519
## siRNAsiPARP1:genotypeALC1KO 0.10889 0.11243 0.969 0.339
## H2O2:siRNAsiPARP1:genotypeALC1KO 0.01371 0.09539 0.144 0.886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1109 on 40 degrees of freedom
## Multiple R-squared: 0.9246, Adjusted R-squared: 0.9114
## F-statistic: 70.09 on 7 and 40 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit3))
## AIC: -65.63346
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ H2O2*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit4))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ H2O2 * siRNA * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 601.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0069 -0.4188 0.1341 0.4384 2.2311
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 3255 57.06
## Residual 127693 357.34
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2054.52 184.07 35.29 11.161 3.97e-13
## H2O2 -1018.50 153.65 32.00 -6.629 1.78e-07
## siRNAsiPARP1 -117.27 260.32 35.29 -0.451 0.6551
## genotypeALC1KO 361.79 260.32 35.29 1.390 0.1733
## H2O2:siRNAsiPARP1 -116.96 217.29 32.00 -0.538 0.5941
## H2O2:genotypeALC1KO -410.10 217.29 32.00 -1.887 0.0682
## siRNAsiPARP1:genotypeALC1KO 423.11 368.15 35.29 1.149 0.2582
## H2O2:siRNAsiPARP1:genotypeALC1KO -91.64 307.30 32.00 -0.298 0.7675
##
## (Intercept) ***
## H2O2 ***
## siRNAsiPARP1
## genotypeALC1KO
## H2O2:siRNAsiPARP1
## H2O2:genotypeALC1KO .
## siRNAsiPARP1:genotypeALC1KO
## H2O2:siRNAsiPARP1:genotypeALC1KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) H2O2 sRNAsPARP1 gALC1K H2O2:sRNAPARP1 H2O2:A sRNAPARP1:
## H2O2 -0.809
## siRNAsPARP1 -0.707 0.572
## gntypALC1KO -0.707 0.572 0.500
## H2O2:sRNAPARP1 0.572 -0.707 -0.809 -0.404
## H2O2:ALC1KO 0.572 -0.707 -0.404 -0.809 0.500
## sRNAPARP1:A 0.500 -0.404 -0.707 -0.707 0.572 0.572
## H2O2:RNAPARP1: -0.404 0.500 0.572 0.572 -0.707 -0.707 -0.809
cat("AIC: ", AIC(fit4))
## AIC: 621.279
simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula
fit5 <- lm(Counts ~ Experiment + poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit5))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 2) * siRNA * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -725.74 -158.61 -17.09 102.58 475.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1021.52 98.18 10.405 4.18e-12
## Experimentexp2 148.69 98.18 1.514 0.1392
## Experimentexp3 -9.75 98.18 -0.099 0.9215
## poly(H2O2, 2)1 -4737.38 555.39 -8.530 5.80e-10
## poly(H2O2, 2)2 61.12 555.39 0.110 0.9130
## siRNAsiPARP1 -230.58 113.37 -2.034 0.0498
## genotypeALC1KO -35.50 113.37 -0.313 0.7561
## poly(H2O2, 2)1:siRNAsiPARP1 -544.03 785.44 -0.693 0.4932
## poly(H2O2, 2)2:siRNAsiPARP1 1647.83 785.44 2.098 0.0434
## poly(H2O2, 2)1:genotypeALC1KO -1907.53 785.44 -2.429 0.0206
## poly(H2O2, 2)2:genotypeALC1KO 1754.59 785.44 2.234 0.0322
## siRNAsiPARP1:genotypeALC1KO 334.33 160.33 2.085 0.0446
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO -426.25 1110.78 -0.384 0.7036
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41 1110.78 -1.525 0.1364
##
## (Intercept) ***
## Experimentexp2
## Experimentexp3
## poly(H2O2, 2)1 ***
## poly(H2O2, 2)2
## siRNAsiPARP1 *
## genotypeALC1KO
## poly(H2O2, 2)1:siRNAsiPARP1
## poly(H2O2, 2)2:siRNAsiPARP1 *
## poly(H2O2, 2)1:genotypeALC1KO *
## poly(H2O2, 2)2:genotypeALC1KO *
## siRNAsiPARP1:genotypeALC1KO *
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.7 on 34 degrees of freedom
## Multiple R-squared: 0.9403, Adjusted R-squared: 0.9175
## F-statistic: 41.2 on 13 and 34 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit5))
## AIC: 689.8121
simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit6))
##
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 2) * siRNA * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.233442 -0.077535 -0.008586 0.088724 0.298066
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.000e+00 4.284e-02 23.343
## poly(H2O2, 2)1 -4.464e+00 2.968e-01 -15.040
## poly(H2O2, 2)2 1.093e-01 2.968e-01 0.368
## siRNAsiPARP1 7.412e-17 6.058e-02 0.000
## genotypeALC1KO -5.228e-17 6.058e-02 0.000
## poly(H2O2, 2)1:siRNAsiPARP1 -1.821e+00 4.197e-01 -4.338
## poly(H2O2, 2)2:siRNAsiPARP1 1.907e+00 4.197e-01 4.544
## poly(H2O2, 2)1:genotypeALC1KO -1.965e+00 4.197e-01 -4.682
## poly(H2O2, 2)2:genotypeALC1KO 1.670e+00 4.197e-01 3.978
## siRNAsiPARP1:genotypeALC1KO 9.821e-17 8.568e-02 0.000
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 1.572e+00 5.936e-01 2.649
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00 5.936e-01 -3.593
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(H2O2, 2)1 < 2e-16 ***
## poly(H2O2, 2)2 0.714841
## siRNAsiPARP1 1.000000
## genotypeALC1KO 1.000000
## poly(H2O2, 2)1:siRNAsiPARP1 0.000111 ***
## poly(H2O2, 2)2:siRNAsiPARP1 6.00e-05 ***
## poly(H2O2, 2)1:genotypeALC1KO 3.95e-05 ***
## poly(H2O2, 2)2:genotypeALC1KO 0.000322 ***
## siRNAsiPARP1:genotypeALC1KO 1.000000
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 0.011915 *
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO 0.000970 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1484 on 36 degrees of freedom
## Multiple R-squared: 0.98, Adjusted R-squared: 0.9738
## F-statistic: 160 on 11 and 36 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit6))
## AIC: -34.74272
simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit7))
##
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 2) * siRNA * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.098930 -0.035263 -0.003283 0.033948 0.125891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.50853 0.01779 28.581 < 2e-16
## poly(H2O2, 2)1 -2.27004 0.12327 -18.415 < 2e-16
## poly(H2O2, 2)2 0.05558 0.12327 0.451 0.654770
## siRNAsiPARP1 -0.12514 0.02516 -4.973 1.63e-05
## genotypeALC1KO -0.12372 0.02516 -4.917 1.94e-05
## poly(H2O2, 2)1:siRNAsiPARP1 -0.13957 0.17433 -0.801 0.428600
## poly(H2O2, 2)2:siRNAsiPARP1 0.71759 0.17433 4.116 0.000214
## poly(H2O2, 2)1:genotypeALC1KO -0.20399 0.17433 -1.170 0.249633
## poly(H2O2, 2)2:genotypeALC1KO 0.62899 0.17433 3.608 0.000930
## siRNAsiPARP1:genotypeALC1KO 0.12217 0.03559 3.433 0.001517
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 0.06375 0.24654 0.259 0.797440
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -0.80896 0.24654 -3.281 0.002301
##
## (Intercept) ***
## poly(H2O2, 2)1 ***
## poly(H2O2, 2)2
## siRNAsiPARP1 ***
## genotypeALC1KO ***
## poly(H2O2, 2)1:siRNAsiPARP1
## poly(H2O2, 2)2:siRNAsiPARP1 ***
## poly(H2O2, 2)1:genotypeALC1KO
## poly(H2O2, 2)2:genotypeALC1KO ***
## siRNAsiPARP1:genotypeALC1KO **
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06164 on 36 degrees of freedom
## Multiple R-squared: 0.9791, Adjusted R-squared: 0.9727
## F-statistic: 153 on 11 and 36 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit7))
## AIC: -119.0964
simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(H2O2,2)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit8))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 2) * siRNA * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 505
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.44832 -0.62194 0.01954 0.46624 1.90599
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 19567 139.9
## Residual 62446 249.9
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1067.83 108.29 8.00 9.861
## poly(H2O2, 2)1 -4737.38 499.78 28.00 -9.479
## poly(H2O2, 2)2 61.12 499.78 28.00 0.122
## siRNAsiPARP1 -230.58 153.14 8.00 -1.506
## genotypeALC1KO -35.50 153.14 8.00 -0.232
## poly(H2O2, 2)1:siRNAsiPARP1 -544.03 706.80 28.00 -0.770
## poly(H2O2, 2)2:siRNAsiPARP1 1647.83 706.80 28.00 2.331
## poly(H2O2, 2)1:genotypeALC1KO -1907.53 706.80 28.00 -2.699
## poly(H2O2, 2)2:genotypeALC1KO 1754.59 706.80 28.00 2.482
## siRNAsiPARP1:genotypeALC1KO 334.33 216.57 8.00 1.544
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO -426.25 999.57 28.00 -0.426
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41 999.57 28.00 -1.695
## Pr(>|t|)
## (Intercept) 9.42e-06 ***
## poly(H2O2, 2)1 3.09e-10 ***
## poly(H2O2, 2)2 0.9035
## siRNAsiPARP1 0.1706
## genotypeALC1KO 0.8225
## poly(H2O2, 2)1:siRNAsiPARP1 0.4479
## poly(H2O2, 2)2:siRNAsiPARP1 0.0272 *
## poly(H2O2, 2)1:genotypeALC1KO 0.0117 *
## poly(H2O2, 2)2:genotypeALC1KO 0.0193 *
## siRNAsiPARP1:genotypeALC1KO 0.1612
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 0.6731
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO 0.1011
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(H2O2,2)1 pl(H2O2,2)2 sRNAsPARP1 gALC1K
## pl(H2O2,2)1 0.000
## pl(H2O2,2)2 0.000 0.000
## siRNAsPARP1 -0.707 0.000 0.000
## gntypALC1KO -0.707 0.000 0.000 0.500
## pl(H2O2,2)1:RNAPARP1 0.000 -0.707 0.000 0.000 0.000
## pl(H2O2,2)2:RNAPARP1 0.000 0.000 -0.707 0.000 0.000
## p(H2O2,2)1:A 0.000 -0.707 0.000 0.000 0.000
## p(H2O2,2)2:A 0.000 0.000 -0.707 0.000 0.000
## sRNAPARP1:A 0.500 0.000 0.000 -0.707 -0.707
## p(H2O2,2)1:RNAPARP1: 0.000 0.500 0.000 0.000 0.000
## p(H2O2,2)2:RNAPARP1: 0.000 0.000 0.500 0.000 0.000
## pl(H2O2,2)1:RNAPARP1 pl(H2O2,2)2:RNAPARP1 p(H2O2,2)1:A
## pl(H2O2,2)1
## pl(H2O2,2)2
## siRNAsPARP1
## gntypALC1KO
## pl(H2O2,2)1:RNAPARP1
## pl(H2O2,2)2:RNAPARP1 0.000
## p(H2O2,2)1:A 0.500 0.000
## p(H2O2,2)2:A 0.000 0.500 0.000
## sRNAPARP1:A 0.000 0.000 0.000
## p(H2O2,2)1:RNAPARP1: -0.707 0.000 -0.707
## p(H2O2,2)2:RNAPARP1: 0.000 -0.707 0.000
## p(H2O2,2)2:A sRNAPARP1: p(H2O2,2)1:RNAPARP1:
## pl(H2O2,2)1
## pl(H2O2,2)2
## siRNAsPARP1
## gntypALC1KO
## pl(H2O2,2)1:RNAPARP1
## pl(H2O2,2)2:RNAPARP1
## p(H2O2,2)1:A
## p(H2O2,2)2:A
## sRNAPARP1:A 0.000
## p(H2O2,2)1:RNAPARP1: 0.000 0.000
## p(H2O2,2)2:RNAPARP1: -0.707 0.000 0.000
cat("AIC: ", AIC(fit8))
## AIC: 533.0262
simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula
fit9 <- lm(Counts ~ Experiment + poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit9))
##
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 3) * siRNA * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -747.69 -119.62 6.48 112.67 451.62
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1021.521 94.245 10.839
## Experimentexp2 148.688 94.245 1.578
## Experimentexp3 -9.750 94.245 -0.103
## poly(H2O2, 3)1 -4737.379 533.132 -8.886
## poly(H2O2, 3)2 61.116 533.132 0.115
## poly(H2O2, 3)3 -743.243 533.132 -1.394
## siRNAsiPARP1 -230.583 108.825 -2.119
## genotypeALC1KO -35.500 108.825 -0.326
## poly(H2O2, 3)1:siRNAsiPARP1 -544.032 753.963 -0.722
## poly(H2O2, 3)2:siRNAsiPARP1 1647.828 753.963 2.186
## poly(H2O2, 3)3:siRNAsiPARP1 9.997 753.963 0.013
## poly(H2O2, 3)1:genotypeALC1KO -1907.534 753.963 -2.530
## poly(H2O2, 3)2:genotypeALC1KO 1754.587 753.963 2.327
## poly(H2O2, 3)3:genotypeALC1KO 106.013 753.963 0.141
## siRNAsiPARP1:genotypeALC1KO 334.333 153.902 2.172
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO -426.247 1066.265 -0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412 1066.265 -1.589
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO -54.356 1066.265 -0.051
## Pr(>|t|)
## (Intercept) 6.78e-12 ***
## Experimentexp2 0.1251
## Experimentexp3 0.9183
## poly(H2O2, 3)1 6.64e-10 ***
## poly(H2O2, 3)2 0.9095
## poly(H2O2, 3)3 0.1735
## siRNAsiPARP1 0.0425 *
## genotypeALC1KO 0.7465
## poly(H2O2, 3)1:siRNAsiPARP1 0.4761
## poly(H2O2, 3)2:siRNAsiPARP1 0.0368 *
## poly(H2O2, 3)3:siRNAsiPARP1 0.9895
## poly(H2O2, 3)1:genotypeALC1KO 0.0169 *
## poly(H2O2, 3)2:genotypeALC1KO 0.0269 *
## poly(H2O2, 3)3:genotypeALC1KO 0.8891
## siRNAsiPARP1:genotypeALC1KO 0.0379 *
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.6922
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 0.1225
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.9597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.6 on 30 degrees of freedom
## Multiple R-squared: 0.9515, Adjusted R-squared: 0.924
## F-statistic: 34.6 on 17 and 30 DF, p-value: 4.022e-15
cat("AIC: ", AIC(fit9))
## AIC: 687.8777
simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit10))
##
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 3) * siRNA * genotype, data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.222920 -0.039175 -0.001219 0.045161 0.187158
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.000e+00 2.783e-02 35.936
## poly(H2O2, 3)1 -4.464e+00 1.928e-01 -23.154
## poly(H2O2, 3)2 1.093e-01 1.928e-01 0.567
## poly(H2O2, 3)3 -6.838e-01 1.928e-01 -3.547
## siRNAsiPARP1 1.666e-16 3.935e-02 0.000
## genotypeALC1KO 1.302e-16 3.935e-02 0.000
## poly(H2O2, 3)1:siRNAsiPARP1 -1.821e+00 2.727e-01 -6.679
## poly(H2O2, 3)2:siRNAsiPARP1 1.907e+00 2.727e-01 6.996
## poly(H2O2, 3)3:siRNAsiPARP1 -1.882e-01 2.727e-01 -0.690
## poly(H2O2, 3)1:genotypeALC1KO -1.965e+00 2.727e-01 -7.208
## poly(H2O2, 3)2:genotypeALC1KO 1.670e+00 2.727e-01 6.124
## poly(H2O2, 3)3:genotypeALC1KO 5.081e-02 2.727e-01 0.186
## siRNAsiPARP1:genotypeALC1KO -3.274e-16 5.565e-02 0.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 1.572e+00 3.856e-01 4.078
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00 3.856e-01 -5.531
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 2.267e-01 3.856e-01 0.588
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(H2O2, 3)1 < 2e-16 ***
## poly(H2O2, 3)2 0.574725
## poly(H2O2, 3)3 0.001227 **
## siRNAsiPARP1 1.000000
## genotypeALC1KO 1.000000
## poly(H2O2, 3)1:siRNAsiPARP1 1.54e-07 ***
## poly(H2O2, 3)2:siRNAsiPARP1 6.30e-08 ***
## poly(H2O2, 3)3:siRNAsiPARP1 0.495028
## poly(H2O2, 3)1:genotypeALC1KO 3.47e-08 ***
## poly(H2O2, 3)2:genotypeALC1KO 7.59e-07 ***
## poly(H2O2, 3)3:genotypeALC1KO 0.853350
## siRNAsiPARP1:genotypeALC1KO 1.000000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.000281 ***
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 4.23e-06 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.560695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0964 on 32 degrees of freedom
## Multiple R-squared: 0.9925, Adjusted R-squared: 0.989
## F-statistic: 281.6 on 15 and 32 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit10))
## AIC: -73.81496
simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit11))
##
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 3) * siRNA * genotype,
## data = dataset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.085467 -0.014959 -0.000465 0.019616 0.071756
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.508535 0.011511 44.179
## poly(H2O2, 3)1 -2.270035 0.079750 -28.464
## poly(H2O2, 3)2 0.055582 0.079750 0.697
## poly(H2O2, 3)3 -0.347722 0.079750 -4.360
## siRNAsiPARP1 -0.125138 0.016279 -7.687
## genotypeALC1KO -0.123719 0.016279 -7.600
## poly(H2O2, 3)1:siRNAsiPARP1 -0.139574 0.112783 -1.238
## poly(H2O2, 3)2:siRNAsiPARP1 0.717592 0.112783 6.363
## poly(H2O2, 3)3:siRNAsiPARP1 0.013414 0.112783 0.119
## poly(H2O2, 3)1:genotypeALC1KO -0.203992 0.112783 -1.809
## poly(H2O2, 3)2:genotypeALC1KO 0.628990 0.112783 5.577
## poly(H2O2, 3)3:genotypeALC1KO 0.104147 0.112783 0.923
## siRNAsiPARP1:genotypeALC1KO 0.122166 0.023022 5.307
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.063748 0.159500 0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -0.808958 0.159500 -5.072
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.003173 0.159500 0.020
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(H2O2, 3)1 < 2e-16 ***
## poly(H2O2, 3)2 0.490862
## poly(H2O2, 3)3 0.000126 ***
## siRNAsiPARP1 9.20e-09 ***
## genotypeALC1KO 1.17e-08 ***
## poly(H2O2, 3)1:siRNAsiPARP1 0.224889
## poly(H2O2, 3)2:siRNAsiPARP1 3.82e-07 ***
## poly(H2O2, 3)3:siRNAsiPARP1 0.906070
## poly(H2O2, 3)1:genotypeALC1KO 0.079901 .
## poly(H2O2, 3)2:genotypeALC1KO 3.71e-06 ***
## poly(H2O2, 3)3:genotypeALC1KO 0.362697
## siRNAsiPARP1:genotypeALC1KO 8.14e-06 ***
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.692050
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 1.61e-05 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.984250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03987 on 32 degrees of freedom
## Multiple R-squared: 0.9922, Adjusted R-squared: 0.9886
## F-statistic: 271.6 on 15 and 32 DF, p-value: < 2.2e-16
cat("AIC: ", AIC(fit11))
## AIC: -158.557
simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(H2O2,3)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit12))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 3) * siRNA * genotype + (1 | UID)
## Data: dataset
##
## REML criterion at convergence: 439.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67826 -0.45216 -0.01205 0.40378 1.92855
##
## Random effects:
## Groups Name Variance Std.Dev.
## UID (Intercept) 22071 148.6
## Residual 52429 229.0
## Number of obs: 48, groups: UID, 12
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1067.833 108.288 8.000
## poly(H2O2, 3)1 -4737.379 457.950 24.000
## poly(H2O2, 3)2 61.116 457.950 24.000
## poly(H2O2, 3)3 -743.243 457.950 24.000
## siRNAsiPARP1 -230.583 153.142 8.000
## genotypeALC1KO -35.500 153.142 8.000
## poly(H2O2, 3)1:siRNAsiPARP1 -544.032 647.639 24.000
## poly(H2O2, 3)2:siRNAsiPARP1 1647.828 647.639 24.000
## poly(H2O2, 3)3:siRNAsiPARP1 9.997 647.639 24.000
## poly(H2O2, 3)1:genotypeALC1KO -1907.534 647.639 24.000
## poly(H2O2, 3)2:genotypeALC1KO 1754.587 647.639 24.000
## poly(H2O2, 3)3:genotypeALC1KO 106.013 647.639 24.000
## siRNAsiPARP1:genotypeALC1KO 334.333 216.575 8.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO -426.247 915.899 24.000
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412 915.899 24.000
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO -54.356 915.899 24.000
## t value Pr(>|t|)
## (Intercept) 9.861 9.42e-06 ***
## poly(H2O2, 3)1 -10.345 2.53e-10 ***
## poly(H2O2, 3)2 0.133 0.89495
## poly(H2O2, 3)3 -1.623 0.11766
## siRNAsiPARP1 -1.506 0.17057
## genotypeALC1KO -0.232 0.82250
## poly(H2O2, 3)1:siRNAsiPARP1 -0.840 0.40919
## poly(H2O2, 3)2:siRNAsiPARP1 2.544 0.01780 *
## poly(H2O2, 3)3:siRNAsiPARP1 0.015 0.98781
## poly(H2O2, 3)1:genotypeALC1KO -2.945 0.00706 **
## poly(H2O2, 3)2:genotypeALC1KO 2.709 0.01225 *
## poly(H2O2, 3)3:genotypeALC1KO 0.164 0.87134
## siRNAsiPARP1:genotypeALC1KO 1.544 0.16123
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO -0.465 0.64585
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1.850 0.07666 .
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO -0.059 0.95317
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("AIC: ", AIC(fit12))
## AIC: 475.7654
simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Final Result
fit <- fit11
output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
rownames(output) <- gsub("siRNA", paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1], sep = " in " )
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1], sep = " " )
# suggested result table
kable(output, row.names = T)
| H2O21 in WT siCtrl |
-2.2700352 |
0.0797498 |
-28.4644590 |
0.0000000 |
| H2O22 in WT siCtrl |
0.0555822 |
0.0797498 |
0.6969573 |
0.4908621 |
| H2O23 in WT siCtrl |
-0.3477220 |
0.0797498 |
-4.3601612 |
0.0001261 |
| H2O21: siCtrl vs. siPARP1 in WT |
-0.1395743 |
0.1127833 |
-1.2375448 |
0.2248890 |
| H2O22: siCtrl vs. siPARP1 in WT |
0.7175920 |
0.1127833 |
6.3625752 |
0.0000004 |
| H2O23: siCtrl vs. siPARP1 in WT |
0.0134139 |
0.1127833 |
0.1189354 |
0.9060700 |
| H2O21: WT vs. ALC1KO in siCtrl |
-0.2039924 |
0.1127833 |
-1.8087113 |
0.0799010 |
| H2O22: WT vs. ALC1KO in siCtrl |
0.6289899 |
0.1127833 |
5.5769793 |
0.0000037 |
| H2O23: WT vs. ALC1KO in siCtrl |
0.1041468 |
0.1127833 |
0.9234246 |
0.3626969 |
| H2O21: siCtrl vs. siPARP1: WT vs. ALC1KO |
0.0637482 |
0.1594996 |
0.3996763 |
0.6920503 |
| H2O22: siCtrl vs. siPARP1: WT vs. ALC1KO |
-0.8089581 |
0.1594996 |
-5.0718495 |
0.0000161 |
| H2O23: siCtrl vs. siPARP1: WT vs. ALC1KO |
0.0031734 |
0.1594996 |
0.0198959 |
0.9842500 |
write.table(output, file = "FigureS3_H2O2_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1KO")
fit <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]
rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype", paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1], sep = " in " )
rownames(output) <- gsub("siRNA", paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1], sep = " in " )
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1], sep = " " )
# suggested result table
kable(output, row.names = T)
| H2O21 in ALC1KO siCtrl |
-2.4740275 |
0.0797498 |
-31.0223630 |
0.0000000 |
| H2O22 in ALC1KO siCtrl |
0.6845721 |
0.0797498 |
8.5839971 |
0.0000000 |
| H2O23 in ALC1KO siCtrl |
-0.2435752 |
0.0797498 |
-3.0542416 |
0.0045206 |
| H2O21: siCtrl vs. siPARP1 in ALC1KO |
-0.0758261 |
0.1127833 |
-0.6723172 |
0.5062076 |
| H2O22: siCtrl vs. siPARP1 in ALC1KO |
-0.0913661 |
0.1127833 |
-0.8101032 |
0.4238645 |
| H2O23: siCtrl vs. siPARP1 in ALC1KO |
0.0165873 |
0.1127833 |
0.1470725 |
0.8839975 |
| H2O21: ALC1KO vs. WT in siCtrl |
0.2039924 |
0.1127833 |
1.8087113 |
0.0799010 |
| H2O22: ALC1KO vs. WT in siCtrl |
-0.6289899 |
0.1127833 |
-5.5769793 |
0.0000037 |
| H2O23: ALC1KO vs. WT in siCtrl |
-0.1041468 |
0.1127833 |
-0.9234246 |
0.3626969 |
| H2O21: siCtrl vs. siPARP1: ALC1KO vs. WT |
-0.0637482 |
0.1594996 |
-0.3996763 |
0.6920503 |
| H2O22: siCtrl vs. siPARP1: ALC1KO vs. WT |
0.8089581 |
0.1594996 |
5.0718495 |
0.0000161 |
| H2O23: siCtrl vs. siPARP1: ALC1KO vs. WT |
-0.0031734 |
0.1594996 |
-0.0198959 |
0.9842500 |
write.table(output, file = "FigureS3_H2O2_Stats_Ref_ALC1KO.txt", quote = F, sep = "\t", row.names = T, col.names = NA)